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' We present the numerically exact ground state energy, efltective mass, and isotope exponents of 

' a one- dimensional lattice polaron, valid for any range of electron-phonon interaction, applying a 

04 , continuous-time Quantum Monte Carlo (QMC) technique in a wide range of coupling strength and 

■ adiabatic ratio. The QMC method is free from any systematic finite-size and finite-time-step errors. 

C^h' We compare our numerically exact results with analytical weak-coupling theory and with the strong- 

I coupling 1/A expansion. We show that the exact results agree well with the canonical Frohlich and 

Holstein-Lang-Firsov theories in the weak and strong coupling limits, respectively, for any range of 
' interaction. We find a strong dependence of the polaron dynamics on the range of interaction. An 

increased range of interaction has a similar effect to an increased (less adiabatic) phonon frequency: 
specifically, a reduction in the effective mass. 

^ ' 

^ , I. INTRODUCTION 
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While qualitative features of polarons wer e well recognized a long time ago and have been described in several 
review papers and textbooks (see Ref . IllO lsl^ for recent publications), there is renewed interest in quantitative studies 
owing; to t he ove rwhelming evidence for polaronic carriers in cuprates, fuUerenes, and manganites (see, for example 

(-H , Under certain conditions the multi-polaron system can be metallic but with polaronic carriers rather than bare 
O ' electrons. There is a qualitative difference between the ordinary metal and the polaronic one. One can account for 
O , the electron-phonon (e-ph) interaction in simple metals by applying Migdal's theoremiS. The theorem shows that 
the contribution of diagrams with "crossing" phonon lines (so called "vertex" corrections) is small if the parameter 
Xhuj/Ep is small, where A is the dimensionless (BCS) e-ph coupling constant, uj is the characteristic phonon frequency, 
^ ■ and Ep is the Fermi energy. Neglecting the vertex corrections, Migdal calculated the renormalized electron mass as 
[ TO* = mo(l -t- A) (near the Fermi level)^^, where mg is the band mass in the absence of electron-phonon interaction, 
. and Eliashbergi^ extended Migdal's theory to describe the BCS superconducting state at intermediate values of A, 
A < 1. Later on many authors applied Migdal-Eliashberg theory with A much larger than 1 (see, for example, Ref. IT^ . 
^— ^ J On the other hand, starting from the infinite coupling limit, A = oo and applying the inverse (1/A) expansion 
■ techniquei^ one can showiSiiLiS that the many-electron system collapses into the small polaron regime at A ~ 1 almost 
independently of the adiabatic ratio fku/Ep. This regime is beyond Migdal-Eliashberg theory, where the effective mass 
"Y^ ' approximation is used and the electron bandwidth is infinite. It is a well established theorem that a self-trapping 
Ci ' crossover is analytical in the coupling strength, so that one could believe that the sum of all diagrams (including the 
vertex corrections) should produce the exact solution if the expansion converges. Indeed, results of QMC simulations 
based on summing the Fcynman diagrams provide the exact answer in the continuous (large polaron) model. On the 
other hand, the small polaron regime cannot be reached by summation of the standard Feynman-Dyson perturbation 
diagrams using a translation-invariant Green function G(r, r',r) = G(r — r',r) with the Fourier transform G(k, il) 
^ [ prior to solving the Dyson equations on a discrete lattice. This assumption excludes the possibility of local violation of 
the translational symmetrySS due to the lattice deformation in any order of the Feynman-Dyson perturbation theory 
similar to the absence of the anomalous (Bogoliubov) averages in any order of perturbation theorjii^. One way to 
describe the formation of the lowest polaronic band is to introduce an infinitesimal translation-noninvariant potential, 
which should be set zero only in the final solution obtained by the summation of Feynman diagrams for the Fourier 
5^ , transform G(k, k', fi) of G(r, r', r) rather than for G(k, il)^^. As in the case of the off-diagonal superconducting order 
parameter, the off-diagonal terms of the Green function, in particular the Umklapp terms with k' = k + G, drive the 
system into a small polaron ground state at sufficiently large coupling. Setting the translation-noninvariant potential 
to zero in the solution of the equations of motion restores the translation symmetry but in a polaron band rather than 
in the bare electron band, which turns out to be an excited stated. Alternatively, one can work with momentum 
eigenstates throughout the whole coupling region, but taking into account the finite electron bandwidth (i. c., including 
Umklapp terms). In recent years many such numerical and analytical studies have confirme d the conclusion^*' that the 
Migdal-Eliashberg theory breaks down at A > 1 (see Refs. l22l23l24l25l26l27l28l29l3Ci3ll32M33ll34M35M36l and references 
therein) . 

In ordinary metals, where the Migdal approximation is believed to be valid, the renormalized effective mass of 
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electrons is independent of the ion mass M because the electron-phonon interaction constant A does not depend on 
M . However, when the e-ph interaction is sufficiently strong, the electrons form polarons dressed by lattice distortions, 
with an effective mass m* = moeyip{"/Ep/huj). Here Ep is the polaron binding energy (or the polaron shift), and 
7 is a numerical constant that depends on the radius of electron-phonon interaction and is typically less than 1. 
While Ep in the above expression does not depend on the ion mass, the phonon frequency does. As a result, there 
is a large isotope effect on the carrier mass in polaronic conductors, am — {l/2)ln{m* /m)^'^ , in contrast to the zero 
isotope effect in ordinary metals. Such an effect was found experimentally in the cuprates^*' and manganites'^^. A 
recent high resolution angle resolved photoemission spectroscopy study'^^ provided further compelling evidence for 
strong e-ph interaction in the cuprates. It revealed a fine phonon structure in the electron self-energy of underdoped 
La2_2;Sr2.Cu04 sample a'^^i"^° and a complicated isotope effect in the electron spectral function of Bi2212 that depended 
on the electron energy and momentumii. 

With increasing phonon frequency the range of validity of the 1/A polaron expansion extends to smaller values 
of Aii. As a result, the region of applicability of the Migdal-Eliashberg approach (even with vertex corrections) 
shrinks to smaller values of the coupling, A < 1, with increasing to. Strong correlations between carriers might 
reduce this region furthcr^^. Carriers in the fascinating novel materials are strongly coupled with high-frequency 
optical phonons, making small polarons and non-adiabatic effects relevant for high-temperature superconductivity and 
colossal magnetoresistance phenomena. Indeed the characteristic phonon energies 0.05-0.2 eV in cuprates, manganites 
and in doped fuUerenes are of the same order as the generally accepted values of the hopping integrals t ~ 0.1-0.3 
eVia. 

The continued interest in polarons extends beyond physical description of low-mobility conductors such as the oxides 
or doped polymers. The field has been a testing ground for analytical and numerical techniques for several decades. In 
the past 25 years, several families of powerful numerical methods have been developed and successfully applied to one-, 
two-, and multi-polaron lattice models. These are the quantum Monte Carlo (QMC) simulationsi2il244iiS4SiiLi2ii2iS, 
exact diagonalization of finite clustersiSM^iiSi^, advanced variational methoda^i^SiSiiMi^, and the density-matrix 
renormalization group^. Many methods have been developed so far as to enable reliable calculation of not only static 
and thermodynamic polaron properties, but also of the effective mass, spectrum, and, in some cases, the spectral 
function of the polaron. 

At the same time, the bulk of the lattice polaron studies have been limited to the short-range electron-phonon 
interactions described by the Holstein modelS. In numerical calculations, the locality of the interaction usually 
simplifies the algorithm and reduces the finite-size errors. However, as pointed out by two of us (AK)'^^, the Holstein 
model is not a typical but an extreme polaron model because the screening length is normally larger than the lattice 
constant in doped insulators. It yields the highest possible value of the polaron mass in the strong coupling limit, if 
lattice vibrations are isotropic or polarised perpendicular to the hopping direction^^. With an on-site electron-phonon 
interaction, during every polaron hop the existing lattice deformation has to relax completely to the undcformed 
state, while a full deformation has to form again at the new location of the particle. Such a process results in the 
exponentially small overlap between the initial and the final states, and in an exponentially large effective mass 
with 7 = 1. Real ionic solids with low density of free carriers are characterized by poor screening and are more 
appropriately described by a long-range electron-phonon interaction. Thus the lattice Frdhlich model introduced in 
Ref. |3a is intermediate between the extremes of the Holstein and FrohlichS limits. On one hand, it is a lattice model 
(like the Holstein one), and the ratio of the hopping integral to the phonon frequency is an important parameter. 
On the other hand, the electron-phonon interaction is long-range, as in the Frohlich model. It was shown in Ref. |3^ 
that in this intermediate case the polaron mass still grows exponentially with the polaron binding energy Ep but the 
parameter 7 is now less than unity. That leads to much reduced numerical values of the polaron mass, hence the term 
mobile small polaron. The model was further studied by numerical cluster diagonalization^ and 1/A expansioniiS. 
In addition, the tw;o-particle model with non-local electron-phonon interactions was studied variationally^^ and by 
the 1/A expansion techniqueii^. These studies confirmed the original conclusion that a long-range interaction 
significantly reduces the effective mass of the carrier, polaron or bipolaron, sometimes by several orders of magnitude, 
in comparison with the Holstein model. It also makes the self-trapping transition more gradual as a function of A, 
and better describable by the Lang-Firsov theory^^. These findings are in agreement with some earlier studies on 
long-range interactions in narrow-band models^*^. 

In this paper, we further generalize the lattice Frohlich model of Ref. [s^ to electron-phonon interaction of some 
finite radius R. We perform a systematic study of the single polaron problem in one dimension as a function of R. In 
the local limit i? — > we recover the results of the Holstein model obtained in the past by various methods mentioned 
above. In the infinite- i? limit the original AK model and its results are fully recovered as well. Our computational tool 
will be the continuous-time path-integral quantum Monte Carlo algorithm developed previously by one of us^. This 
method is particularly suited for investigating long-range electron-phonon interactions because the phonon degrees of 
freedom are integrated out analytically. Thus the shape of the interaction does not complicate the algorithm at all, 
but simply modifies the weight function of a Monte Carlo configuration. The method works on infinite lattices and in 
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arbitrary dimensions, eliminating finite-size errors; there is also no truncation of the phonon Fock space. The method 
is also free from finite-time-step errors because it is formulated in continuous time. The method enables unbiased 
calculation (i.e. no numerical errors besides statistical fluctuations) of the polaron energies, effective mass, spectrum, 
density of states, isotope exponents, the number of excited phonons, and other quantities. 

In addition to presenting novel results on the finite-radius Frohlich model, we use the present paper to explain 
many technical details of the polaron QMC method^, which have not previously been published. The electron- 
phonon Hamiltonian is introduced in section II. In section III we describe the continuous-time Monte Carlo method. 
In sections IV- VI we present the numerical results for the energy, effective mass, number of dressing phonons, and 
isotope exponents of lattice polarons for different R, and compare them with weak-coupling and strong-coupling 
analytical results and with numerical results of other authors. Section VII summarizes our conclusions. 



II. ELECTRON-PHONON MODEL 



A. General model Hamiltonian 



The electron-phonon model under investigation represents a single electron interacting with all the ions of an infinite 
hypercubic lattice, with one vibrational degree of freedom per unit cell. The Hamiltonian takes the form 

H = Hc + Hph + i^c-ph, (1) 

where 

iJe=-t^4c„,, (2) 
(nn'> 
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m m 

and 

-ffo-ph - X! /m(n)c],c„^m- (4) 

mil 

The free-electron term He describes the movement of a single electron through the lattice by the process of nearest- 
neighbor hopping. Here the operator cj, creates an electron on site n, the operator Cn' destroys an electron on site 
n', and (nn') denotes pairs of nearest-neighbor sites. The phonon term Hph represents the vibrations of the lattice 
ions. Here the operator is the displacement of the mth ion from its equilibrium position, and = —ihd/dS,m its 
momentum. It is assumed that the ions, each of ionic mass Af , are non-interacting and so have the same characteristic 
(phonon) frequency lo. The final part of the Hamiltonian, the electron-phonon term Hc^ph, is of the "density- 
displacement" type, where the interaction energy between the electron and the mth ion is proportional to (the 
displacement of the mth ion from its equilibrium position). Here cJ^Cn is the electron number operator, and /m(n) is 
interpreted as the interaction force between the electron on site n and the mth lattice ion. 

The model is parameterized by two dimensionless quantities. The first is the dimensionless phonon frequency 

Co — huj/t. (5) 

The second is related to the small-polaron binding energy Ep, derived in section (|Vip . which serves as a natural and 
convenient measure of the strength of the electron-phonon interaction. The dimensionless electron-phonon coupling 
constant is defined as 



where zt is the bare-electron half bandwidth, with z the lattice coordination number. 
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FIG. 1: Geometry of the Frohlich model l[THHl shown in one dimension. The mobile charge carrier moves on the lower chain 
with nearest-neighbor hopping integral t and interacts with all the ions of the upper chain. The displacements of the ions 
are polarized in a direction perpendicular to the chains. 



B. Discrete Frohlich model 



Some time ago Alexandrov and Kornilovitch proposed a long-range discrete Frohlich interactionp^ to describe the 
interaction between a hole and the apical oxygen ions in high-Tc superconducting materials. The model is depicted in 
FIG. ^for the one-dimensional case. The mobile carrier (electron or hole) may hop from site to nearest-neighbor site 
along the lower chain. The chain consists of an infinite number of lattice sites with lattice constant a. The electron 
interacts with all the ions which reside at the lattice sites of a similar chain that is parallel to the first. The separation 
of the two chains is equal to the lattice constant a. We assume that the vibrations of the ions are polarized in a 
direction that is perpendicular to the chains, and that the ions do not interact with each other. 

Let us find the appropriate form for the interaction force fmin) between the mobile charge-carrier on the nth site 
(of the lower chain) and the mth ion (of the upper chain). Since both m and n are measured in units of a, we choose 
from this point on to take a — 1. The presence of the charge-carrier displaces the mth ion by a small distance in 
a direction perpendicular to the chain, as shown in FIG.^ By expanding the Coulomb potential in powers of ^m, we 
deduce2& that the Hamiltonian for the discrete Frohlich model is that of our generalized model Hamiltonian, Eq. ^ 
with the electron-phonon interaction force having the form 

fm{n) = ^ (7) 

with a constant k. Physically, this model was proposed in order to represent the interaction between a hole in 
the copper-oxygen layer (lower chain) and the apical oxygens in the ionic layer (upper chain) contained within the 
structure of certain doped high-Tc superconductors such as YBa2Cu30g+2i^. These materials are highly anisotropic 
due to the fact that the holes are sharply localized in the copper-oxygen layer, giving rise to poor conduction in the 
c-direction (normal to copper-oxygen layer). This leads to very poor screening of the electron-phonon interaction 
in the c-direction, and almost complete screening in the a-h plane. This justifies the restriction to phonon modes 
polarized in the c direction. 



C. Screened Frohlich model 



Our aim in this paper is to investigate the way in which the shape of the long-range electron-phonon interaction 
affects the properties of the polaron. It is therefore interesting to study the screened Frohlich model, in which 
the screening effect due to the presence of other electrons in the lattice is taken into account from within fm{n). 
Accordingly, let us define the interaction force for the screened Frohlich model as 

where R^c is the screening length. That is, the screened force is the unscreened force multiplied by an exponential 
damping factor. Increasing the value of R^c decreases the screening effect and thus increases the width of the interaction 
force. 
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FIG. 2: The shape of the screened Frohhch interaction force, Eq. |H| at screening lengths of R^c (bold line, Holstein 
interaction), i?sc = 1 (dashed line), Rac = 3 (dot-dashed line) and iiac oo (thin line, non-screened Frohlich interaction). 



The Holstein model describes an electron that interacts only with the oscillator it currently occupies ( "short-range" 
interaction). This may be regarded as a special case of equation ^ with Rsc — > 0, so that fm{n) = —KSmn- Simply 
by altering the value of the parameter Rgc, we can easily cross over from the Holstein model, through the screened 
Frohlich model, to the unscreened Frohlich model, in a universal manner. In this paper we consider the following four 
cases in one dimension: 

1. Holstein model with Rsc 0. 

2. Screened Frohlich model, with Rsc = 1- 

3. Screened Frohlich model, with Rsc = 3. 

4. Unscreened Frohlich model with Rsc oo. 

The shapes of the electron-phonon interaction force fmin) for each of the above screening lengths are shown in 
FIG. 13 Note that, based on calculations involving the dynamic properties of the polaron response-^^, the amount of 
screening we impose here is greater than that expected in the high-Tc compounds. 



III. PATH INTEGRAL APPROACH 



A. Effective mass using a partial partition function 

The effective mass of the polaron m* is defined for the isotropic or one-dimensional case as^ 

where E(P) is the ground state energy for total momentum P (sum of the momenta of the electron and all the 
phonons). 

The evaluation of (jn*)^^ by differentiating QMC energies is not practical within our approach because a minus- 
sign problem arises for finite momentum, exacerbating the errors already present in such a procedure. The usual 
means of extracting dynamical properties (such as the effective mass) from QMC simulation is by making use of some 
kind of analytical continuation from imaginary to real time. However, it is possible to infer m* directly from QMC 
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simulation by considering electron trajectories with twisted (rather than periodic) boundary conditions in imaginary 
time. Kornilovitch'*^ showed (for the isotropic or one-dimensional case) that 

where toq = / {2ta?) is the bare electron mass, and 

^Ar = j d^C({U+r'-r}, v' \ e~^" \ {^j, r) , (11) 

is a "partial partition function" (which is similar in form to the total partition function of the system), with 



(12) 



Here |r) is the electron basis, |{fm}) — 1^17^21^3: ' ' ' i^iv) is the ionic displacement basis, and the summations over 
Ar includes all possible values of Ar = r' — r. 

Given Eq. the effective mass may be obtained from QMC simulation by taking the statistical average of (Ar)^, 
sampled over trajectories of the path integral formed from ^Ar- The dissimilarity between the "bra" and "ket" states 
in Eg. II II produces a path integral of ^Ar having twisted (rather than the usual periodic) boundary conditions. Note 
that, since we need only consider the case of P = 0, there is no sign problem. 

B. Continuous imaginary time 

QMC schemes have recently been developed that are implemented directly in continuous imaginary time for lat- 
tice modelsi^iSMij eliminating the problematic finite-time-step error associated with the traditional discrete-time 
approach. The partial partition function Zat in Ea. llll is given in continuous- imaginary-time path- integral form as^ 

^Ar- / V^Vrexp{S), (13) 

J tw 

where the phonon action reads 



^ ^ "'0 
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(14) 



with S,in{T) — dS^m{T)/dT. In forming the path integral above, an imaginary-time dimension r has been introduced, 
having the range < r < The electron and phonon coordinates are represented as continuous functions of 
imaginary time, r(r) and ^m(''"), which can be interpreted as continuous trajectories in t. 

The symbol J^.^ in Ea. ll3l represents the integration over all possible trajectories under twisted boundary conditions 
in imaginary time. The "end states" of the individual trajectories, which arc identified with the states ({Cm+r'-r}i r' | 
and I {fm}, r) in Eq. ^2 are given by^ 

|{^„(0)},r(0)) = 

I Um (/?)}, r(/3)) - I Um+Ar},r + Ar), (15) 

that is, the final state (r = (3) is the initial state (r = 0) with all the coordinates (electron and all phonons) shifted 
by Ar. 

We may decompose the trajectory ^m(T) in Eq.^jinto the sum of a the classical path (the trajectory that extremizes 
5m) and a deviation (or "quantum fluctuation") from it. The part of S'm that contains no terms involving quantum 
fluctuation is the classical action. The classical action is an important quantity, and is given by^SiSS 

St = 2ns!!Mnu:p) {-Km(0) + d(/3)]cosh(fi^/3) + 2em(0)U(/3)} +^m(0)Sm(r) +^m(/3)Cm(r) 

+ ^ fj^dr' fMr))G{T, r')fMr')), (16) 
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where 



c„(.) . (18) 



Q sinh(?ki;/?)' 



and the Green function is 



^, i\_ \ J sinh(fia;r) sinh[fta;(/? — r')], < r < r' , ^ 

^ ' ~ kujsin\,{nw[3) \ sinh[fi^(/3 - r)] sinh(nt^T'), r' < r < /? ' ^^"^^ 

Note that the phonon coordinates in Eq. ^]are those of the end-points only: Cm(0) and Cm(/3)- 



C. Analytical Phonon Integration 

We wish to integrate out the phonon degrees of freedom from the problem analytically, that is, perform the phonon 
path integral 



[ P|exp(^5„,) =ctw / d^exp(^5^). 



(20) 



where the non-classical part of S (terms involving quantum fluctuation) integrates to an unimportant constant^ Ctw, 
reducing the problem to the integration of the classical action S'^. The integration must be performed under twisted 
boundary conditions in imaginary time. Accordingly, we impose the constraints 

em(0)-U ; Cm(/3)=U-Ar (21) 

on S'^^in Eg. [TCI which one can see produces mixed variable terms involving ^mCm-Ar- The phonon integration cannot 
directly be performed in this form. However, we may proceed by transforming into real Fourier components Oq 
and 6q: 

U = ^EK + *Me'''"' (22) 

where D is the dimensionality of the lattice and DN is the total number of phonon degrees of freedom. In this 
representation, the transformed action S'^ is diagonal, and so the phonon path integral in Eg . 1201 decomposes to the 
product of single variable integrals according to 

/tw = Ctw n / da^dh^ exp {S"^) . (23) 

q 

After performing the Gaussian integration (in Cq and 6q) the result is 

7r/isinh(fta;/3) -1-0^/2 



where 



Mw[cosh(fiw/3) - cos(q • Ar)] 
;isinh(nw/3)^^S„(C™_Ar-C„,) , nsinh(fii^/3)^„(B„ + C, 



exp A, (24) 



2Mw[cosh(fi,cj/3) - cos(q • Ar)] AM uj[cos\i{huj (3) - cos(q • Ar)] 

+ E ^ f fdTdT'f^{v{T))G{T, T')f{v{T')), (25) 



which does not contain any phonon degrees of freedom. We have thus transformed the problem from that of an 
electron interacting with many phonons to that of an electron with retarded self-interaction, which allows the QMC 
method to be applied effectively. 
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D. Low temperature limit 

The result above may be conveniently rendered into the required low temperature limit using cosh hw(3 « sinh hjjj(3 « 
ifjiup ^ 1 to give 

/tw cx exp A = exp (Apcr + AA) , (26) 

where 

is the low temperature action for periodic boundary conditions, and 

= j'j'^dfdf'e-^^e-^^P-^'^ [$Ar (r(f), r(f')) - <i>o (r(f), r(f'))] (28) 

is the correction for twisted boundary conditions, in dimensionless form. Here f = It and (3 = t(3 define dimensionless 
imaginary time; the parameters w and A are defined in Eqs.[Sland|n|respectively; and the lattice summation is defined 
as 

$Ar(r(f),r(f')) =^An(r(f))/„,+Ar(r(f')). (29) 

rii 

Note that the dimensionless quantity /m(n) represents the shape or form of the electron-phonon interaction force, 
defined via the decomposition 

/m(n) = kJM. (30) 

_ 1/2 

where k = \^z\Mt'^uj^ / (^^ X^m /m(0))] takes the dimensions of force. 



IV. CONTINUOUS-TIME MONTE CARLO 



A. Algorithm 

Traditionally, path-integral QMC simulation is implemented in discrete imaginary-time, where the trajectory is 
represented by the position of the electron in each of a large number of imaginary-time slices. The use of discrete time 
introduces the problematic finite-time-step systematic error, which scales with the square of the time-slice- width. 

A path-integral QMC scheme implemented directly in continuous imaginary time has been developed for systems 
with a discrete basisiSiSi. Here, the electron trajectory is represented as finite intervals of imaginary time in which 
the system remains in a particular state, separated by sporadic transitions from one state to another (an electron 
hop). The points in imaginary time at which the state of the system changes are called "kinks", as shown in FIG. 13 
It is necessary to consider the statistics governing different directions of kink independently of one another. For our 
one-dimensional case with nearest-neighbor hopping, we need only consider single left and right kinks. The use of 
continuous time completely eliminates the finite-time-step error, rendering the scheme "numerically exact" . 

If Ng is the number of kinks of direction s, we wish to generate random states according to the Monte Carlo weight 



w{{Ns}) = w,x[{Ns}) «;ph({A^.}), (31) 
where the weight from the electron subsystem 

m'^^e-'p 



is given by the Poisson distribution, and the phonon-induced weight 



Wpi,{Ns) - cxp[A(7V,)] 



(33) 
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FIG. 3: Illustration of a one-dimensional electron trajectory in imaginary time. The point in imaginary time at which the 
electron hops to a neighboring site is known as a kink, (a) Shows a trajectory with three kinks: occurring at n, T2, and T4. 

(b) The same trajectory, but with a kink added at time rs. The entire trajectory above ra is shifted by one lattice parameter. 

(c) The same trajectory as in (a) but with the kink at T2 removed. Again, the entire trajectory above T2 is shifted accordingly. 



is given by Eq. [SEI 

Proposed changes to the shape of the trajectory are generated by the addition or deletion of single kinks. This is 
sufhcient in practice. In order to increase efficiency one might also consider changing a kink-direction, repositioning a 
kink in imaginary-time, or altering the temporal ordering of the kinks. The Metropolis method^iS^ accepts or rejects 
the trial change from state fj, to state fi' with a transition probability P(/i — > = g{fi — > /i')a(/i ~* fi'), where 
— > fj,') is the sampling distribution and a{fi /i') is the acceptance probability. For the case of iV,, > 1 (one or 
more kinks exist), we choose g{Ns -I- 1 ^ Ng) = g{Ns + 1) = 1/2, and so the acceptance probability is given by 

aMN.^N. + l) . min(l,4^^±l^^^W+l)j (34) 
= muJl,j;^exp[A{Ns + l)-A{Ns)]\ (35) 



to add a kink, and similarly 



,{N, + l^Ns) = min <j 1, ^^^^ exp[AiNs) - A{Ns + 1)] |> (36) 



to remove a kink of direction s. For the case of = 0, we can only add a kink, and so g{0 ^ 1) = 1, which gives 

aadd(0 ^ 1) - min |l, ^ exp[^(l) - ^(0)]| (37) 



and 



arcm(l ^ 0) = min |l, ^ exp[A(0) - A(l)]| . (38) 
The continuous-imaginary-time QMC step used in this work has the following structure: 

1. Randomly select a kink direction s for the trial change. In the case of a one-dimensional system, this is left or 
right. 

2. Propose a change to the trajectory by randomly selecting whether to add a new kink (at a random imaginary 
time) or to remove an existing kink (selected in a random fashion). This is done according to the selection 
probabilities g{Ns + 1 ^ N^) and g{Ns Ns + 1). 
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3. Accept or reject the proposed change with probabiUty aadd(^ — > ^ + 1) if adding, or ai-cm{Ns + 1 A^) if 
removing a kink. 

4. If the change has been accepted, then the entire trajectory that hes "above" the kink (i.e. from the imaginary 
time of the kink to /?) is shifted across accordingly. If the proposed change has been rejected, then the trajectory 
is left untouched. 

B. Analytical integration over kinks 

The Metropolis algorithm requires the action, which involves a double integration in imaginary time, to be computed 
on each Monte Carlo step. The fact that the trajectory consists of a series of single kinks, between which the trajectory 
is a straight line (r(f) is constant), leads us to decompose the the action A in Eq. |2Hlinto segments according to 

N,+ l 

A = £%+2g £ A,,, (39) 
j=i j=i k=j+i 

where j and k label the kinks (along trajectories corresponding with r and r' respectively in the double integration), 
such that fj is the imaginary-time at which the j'th kink occurs, with tq = and tn^+i = P- We treat the diagonal Ajj 
and off-diagonal segments Ajk separately in order to increase efficiency. Each Ajk involves the range of imaginary-time 
between successive kinks of Tj_i < f < fj and fk-i < t' < fk, in which the electron coordinate is fixed at r(-f ) = r(Tj_i) 
and r(-f') = r(ffc_i) respectively. Thus the value of the lattice summation <I>Ar(r('rj-i), r(-ffc_i)) given by Eg. 1291 has a 
constant value, allowing the double integration appearing in Eas. l27l and l28l to be treated analytically for each segment 
Ajk- The result after rearrangement reads 



J2 [4^)<i>("'+AW$^r'-)] + E E [A2<^^^'"^ + A,<^t"^] , (40) 



^0 k J = l j=l k=] + l 

where we have used the shorthand $2r^^ ~ ^^'^ i"(^fe-i)) foi' lattice summation defined in Eq. 123 and 

1 r 



and 



where we have defined 



^{tj - fj-i) - K^'^l , (41) 

A2 = iifW^We-^^^^-i-^^"), (43) 

UJ 

A3 = iifO')i^«e-'^(^~^'=+^^-i), (44) 



ifO") = 1 _e-i'(f.-^.-i). (45) 

The action can thus be computed efficiently using this double summation over kinks. 

For the models studied in this paper, /m(n) depends only on the relative lattice distance |m — n|, and tends to 
zero (or is zero) at large distance. Consequently, the lattice summation, Eq. 1291 is a function of the single variable 
r' = r2 — ri — Ar only, namely 

$(r') = E/(m')/(m'-r'), (46) 



which can be evaluated for all possible values of r' in advance of the simulation proper, improving the efficiency of 
the QMC scheme. 
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C. Physical Observables 



We consider four physical observables: the ground state energy, the number of phonons in the polaron cloud, the 
effective mass, and the isotope exponent on the effective mass. For a given observable Q, the expectation value is the 
statistical average over trajectories at P = (ground state), which can be written in the form 



(Q>c 



(47) 



where the phonon degrees of freedom have been integrated out, and w{Ns) is given by Eg. 1311 This corresponds to a 
simple arithmetic average within the QMC simulation. 
The ground state (P = 0) energy estimator is given hy^ 



(48) 



which follows from the corresponding finite-imaginary-time energy estimator**^ in the continuum limit. Within the 
QMC simulation, then, we must gather separate statistics for the total number of kinks J2s ^^"^ ^^'^ quantity 
OA/ dp. One can see that the expression for dA/djS is easily obtained by analytically differentiating Eg . 1401 with 
respect to (5. 

The number of phonons in the polaron cloud iVph quantifies the amount of lattice deformation associated with the 
polaron. The value of A^ph is given by the expectation value of the phonon number operator, which can be isolated 
from the model Hamiltonian using the fact that dH/d{hLu)\^^ = Em'^Lc'm- This can be related to the action via 
the free energy Fq — ~(3^^ InZo to give 



ph 



7 , dlj^dm 



dFo 



d{hLu) 



1/ OA 



(49) 



where dA/diI)\^^ is easily obtained by differentiating Eq. 001 with respect to uj holding the product \uj constant. 

As discussed in section (IIII All , by imposing twisted boundary conditions in imaginary time, dynamical properties can 
be inferred directly from QMC simulation. The effective mass of the polaron m* , for the isotropic or one-dimensional 
case, may be measured using 



Too 

TO* 



21, X'^-' 



(50) 



where the difference in position of the endpoints of the trajectory Ar — r(/3) — r(0) is measured in units of the lattice 
constant a, and toq is the bare electron mass. 

The isotope effect is most often observed via its influence on the superconducting transition temperature T^. The 
dependence of Tc on the mass of the lattice ions M has been found empirically to be Tc oc where a is 

known as the isotope exponent on Tc. In a similar way, let us define the isotope exponent on the effective mass am' 
as 



M dm* 
TO* dM 



-M 



TO* d / Too \ 
Too dM VTO* / 



(51) 



for the isotropic or one-dimensional case. On substitution of the derivative of too /to*, as in Eg. 1501 with respect to 
M, we have 



1 



(Ar) 



(Ar)^ 



:dA\ _ 



{Ary 



dAV 



(52) 



where the ionic mass enters our formalism via the phonon frequency w = {k/M)^/'^ , where k is some "spring 
constant", and dA/du is easily obtained from Eq. ^Uj (For a general Z?-dimensional system, we may define the 



d'th component of the isotope exponent on the effective mass as a 



(d) 



moM/m*i = (2/3)~ 
replaced by Ar^). 



(Ar^i) ) , with Ard — rd{P) ~ r^iff), and thus we may write a^; as Eq. with every Ar 



-M{ml/rnQ^d)d/dM {niQ^d/ml), where 
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D. Simulation details 

The QMC scheme is based on the shxiulation of a single trajectory r(-f ) in imaginary time < f < p. The standard 
MetropoUs algorithm is used to alter the shape of the trajectory by the addition and deletion of single kinks, as 
described in section IjlV Ap . The start of the trajectory r(0) does not change throughout the simulation, but the 
other end r(/3) is "free" (open boundary conditions in imaginary time). In practice, the shape of the trajectory 
was represented using a list containing the imaginary-time, and the direction, of each kink. In addition, we also 
kept track of the value of Ar = r(/3) — r(0), and the the total number of kinks of each direction {A^s}. The major 
computational task is the evaluation, on each Monte Carlo step, of the action given by Eq. (|40|l . (The number of 

exponential-function-evaluations was reduced by storing the values of a[^^ and K^^'' along with each kink, reducing the 
overall computational effort). In order to calculate the expectation values of the observables given in section (|IV Cp . 
separate statistics for the quantities {^gNs)Q, {dA/df})o, {dA/dui\.^-)o, ((Ar) )o, {dA/dui)o, and ((Ar) dA/du})^ 
were gathered every 10-50 Monte Carlo steps. We only consider the case of P = 0, corresponding to the ground state 
of the system, where there is no sign problem. 

The four one-dimensional interaction models studied differ only in the value of the screening length Rsc- The model 
dependency enters the simulation via <I>(r') given by Ea. 1461 For each model, simulations were conducted for various 
different values of the dimcnsionless parameters Q and A. 

The value of (3 was set at a sufficiently large value to enforce the low temperature limit exp((I;/3) oo. For the 
present simulations, a value of /3 > 15 was found to make the finite-temperature error negligible. (Reducing the value 
of u), or increasing A, beyond those studied here would require this value of /3 to be increased). Increasing the value 
of P increases the "length" of the trajectory (which involves more kinks), and thus increases the computational effort 
required to perform each Monte Carlo step. 

For each set of model parameters (-Rsc, and A), between 3 and 6 statistically independent Monte Carlo runs were 
performed, each using a different value of the inverse temperature 15 < /3 < 25. The number of Monte Carlo steps in 
each run was taken to be about five times the "warm-up" period. Typically, the runs consisted of between 1 x 10^ 
and 5 x lO'' steps in total. (The statistics gathered for each set of runs were viewed together graphically, in order to 
better estimate the point at which equilibrium had been reached). Only those statistics gathered after the estimated 
warm-up period were included in the averages. Given that the finite-temperature error is small, and in the absence 
of systematic finite-size and finite-time-step errors, the main source of error is statistical. The size of the statistical 
error depends on /m(n) and (3, as well as on Co and A. For each set of model parameters (i?sc, ^ and A), we performed 
a sufficient number of runs to ensure that the Monte Carlo averages were determined to a statistical error of less than 
one percent. 

V. LIMITING CASES 

A. Strong Coupling (Small Polaron) Regime 

When the electron-phonon coupling is strong, the electron becomes "trapped" in a potential well created by the 
induced lattice distortion. In this case the "size" of the polaron state can become comparable with the lattice constant, 
and the term "small polaron" is used. The condition for small polaron formation is A = Ep/zt > 1, which is referred 
to as the strong-coupling regim^S,. The small polaron can move from site to site (at zero temperature) through the 
action of zero-point motion. 

An analytical method to determine the effective mass and energy dispersion in the strong-coupling regime, for 
a lattice polaron with a general electron-phonon interaction force fmi^)^^, is based on the Lang-Firsov canonical 
transformationi^ (which renders the transformed Hamiltonian diagonal for A oo), followed by a second-order 
perturbation technique that uses 1/A as a small parameter—. For nearest-neighbor hopping only, and forms of fm{n) 
that depend on the relative lattice distance |m — n|, the result for the lowest energy levels reads'^^ 



k',{nq} 

where the summation is over intermediate states with one or more phonons, X]q"q polaronic energy-level 

shift 



(k, I J2nn' *nn' exp 



(/„,(n)-/„(n'))(e"'-"4-e- 



k',W) 



(53) 



(54) 
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(used to define A) corresponds with the sohition for A ^ oo, and the small-polaron dispersion 

£p(k)=t^e-«'We-'''-, (55) 

where is known as the (zero temperature) mass renormalization exponent, given by 

ff'(l) - E [/-(O) - /-(0)/m(l)] • (56) 

m 

The final term of Eq. 1531 is quadratic with respect to the bare hopping amplitude inn' , and produces a negative 
correction to the energy. (Here, riq is the number of phonons with wave vector q, and | k',nq) is an excited state 
that consists of a single electron with wave vector k' and one or more phonons.) It is of order of 1/A^ and is almost 
independent of k2^. 

In the ground state of the system (k = 0), the value of ep(k) is real and small compared with — i?p, and so Eg. 1531 
is dominated by —Ep. Thus, the (dimensionless) strong-coupling result for the ground state energy is 

^ = .A, (57) 
and since — N-piJuo, it follows that the number of phonons in the polaron cloud at k = is given by 

A^ph = ^. (58) 

Assuming that the third term in Ea. l53l is completely independent of k, then the inverse effective mass, for the isotropic 
or one-dimensional case, is given by 



Too d^E{k) 
- mo 



TO* " dk^ 
This can be conveniently expressed in dimensionless form as 



-e-»«. (59) 

A:-»0 



— = exp , (60) 



TO, 

where we have defined the constant 

_ -, _ J2m /m(0)/m(l) /„ x 

which depends only on the shape of the electron-phonon interaction force. The isotope exponent on the effective mass 
am*, defined in Eg. 1511 can be written in terms of w cx M~^/^ as 



TO* LU d /TOqX 

Too 2 du) \ TO* / 



(62) 



for the isotropic (or one-dimensional) case. Thus, the (dimensionless) strong-coupling isotope exponent is given bji^i 

a„i' = (63) 

The value of the dimensionless constant < 7 < 1 must in general be determined numerically: the value calculated 
for each of the one-dimensional interaction models studied in this work is presented in TABLE ^ 



B. Weak Coupling (Large Polaron) Regime 



In 1950 Frohlich et al^ considered the ground state of a polaron in the weak-coupling limit A ^ 1, using second- 
order perturbation theory, under the condition that there is never more than one phonon virtually excited, and the 
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Interaction Model 


7 


Holstein 




1.000 


Screened Frohlich {Rsc ~ 


1) 


0.745 


Screened Frohlich (-Rsc ~ 


3) 


0.531 


Non-screened Frohlich 




0.387 



TABLE I: Values of the dimensionless parameter j — 1~ /in(0)/m(l) / X^m' /m' (0) the one-dimensional models studied 
in this work, where 7 depends only on the shape of the interaction force. 



discreteness of the lattice is unimportant (because the band-electron-like "large polaron" stated is much larger than 
the lattice constant). In our case we use the tight-binding dispersion 

£o(k) = -2tcos(fca) (64) 

rather than the parabolic approximation £o(k) — h^h^ /2mo + 0{'k*), where mo — /{2ta?) used by FrohlichSl. Thus 
we write 

// = eo(k)+;iw5I4^q + ^<=-ph' (65) 
q 

where the electron-phonon term iJe-ph is a small perturbation. Assuming that /m(n) depends only on the relative 
lattice distance |m — n|, then i?e_ph given by Eq. 2|may be written in momentum representation as 

Hc-ph ^-^^h (c[-qCkrflj + h-c.) , (66) 



where we have defined 



and 



(67) 



E™/i(0)' 



r 

Here, the summation is over all values of r = m — n for which /m(n) = /(m ~ n) = /r(0) is a function of a single 
variable only. Using standard second order perturbation theory, with an initial state that consists of an electron of 
momentum fik and no phonons, and an intermediate state that consists of an electron with momentum ?i(k — q) and 
a single phonon of momentum Sq, the energy measured from the bottom of the electron band is given by 

£;(k)=eo(k)-^5]^, (69) 
q 

where we have defined 

H/ = eo(k-q) + ;iw-£o(k). (70) 

Here the ground state energy Eq occurs at k = and is a suitable Brillouin zone average. 

The number of virtual phonons in the polaron cloud is defined as iVph = X]q(0' I '^q^q I 0') where | 0') represents 
the eigenstate for the perturbed Hamiltonian. Using standard first order perturbation theory, this is given by 

^P^ = «E^- (71) 



q 



The effective mass is easily found by differentiating Eq. EH twice with respect to k, according to Eq.jSHI to give 



^ ^ 1 ^ uq I y'^"^ - j .72^ 

m* 9„/2 W3 ' I 



/q 1 2 (2W^'2 - W"W) 



q 
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where W = dW/dh and W" = 9^W/9k^. Differentiating this expression with respect to w, according to equation 
(|(j2|l . gives the isotope exponent on the effective mass as 



4ai2 mo 



E 



q L 



fq\^ W"{W -2huj) I /q |2 2VK'2(W^-3fia;) 



M/4 



The above expressions for the observables may be written in the form 

-2-ATb„(/», 



1 



l-AT„.(/"cS), 



mo 

m* 



and 



"m* = ^— ^Ta^. (/, w), 

mo/m* 



(73) 

(74) 
(75) 
(76) 

(77) 



where T{ f,uj) is the coefficient of the linear term in A, which may be calculated numerically. We are interested in the 
ground state properties and thus must evaluate the above expressions at k ^ 0. The ground state values of T(/, w) 
for each one-dimensional interaction model are presented in TABLE ITU for uj = 1. 



Interaction Model 


E^/m(o) 










Holstein {R^c 0) 




1.000 


0.894 


0.537 


0.537 


0.125 


Screened Frohlich {Rsc = 


= 1) 


1.034 


1.080 


0.736 


0.625 


0.200 


Screened Frohlich {Rbc = 




1.133 


1.257 


0.938 


0.678 


0.266 


Unscreened Frohlich {Rsc - 


oo) 


1.269 


1.394 


1.104 


0.687 


0.306 



TABLE II: The numerically calculated values of the coefficient of the linear term in A in Egs. 1741)771 for each one-dimensional 
interaction model at a) = 1 in the weak coupling limit. For example, the ground state energy H74|l of the Holstein polaron is 
given by Eo/t = -2 - 0.894A. 

The coefficients can be expressed in closed form for the simplest case of the one-dimensional Holstein interaction. 
Here we find that the energy is 



so that 



£;(k) = -2tcos{ka) 



2u)t\ 



2_^ll/2' 



[{uj + 2cos{ka)Y -A] 



2uj 



and 



TArph(/,w) 



T,„.(/,w) 



Ta„.(/"<:^) 



' (w2-f 4w)'/'' 
2w(<I> + 2) 

~ (w2 + 4^)3/2 ' 

2w(a; + 2) 

~ (^2 + mfl'' ' 



(w2 4w) 



5/2 



(78) 



(79) 



(80) 



(81) 



(82) 
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FIG. 4: The variation of ground state energy Eo{0) (triangles), together with potential energy (diamonds) and kinetic energy 
(squares), with coupling A for the one- dimensional Holstein model with a dimensionless phonon frequency of ui — 1. The dashed 
line is the strong coupling perturbation (SCP) result H57|l and the dotted line the weak-coupling perturbation (WCP) result 
(I74I I. One can clearly distinguish the large-polaron, transition and small-polaron regions. 

VI. RESULTS 
A. Holstein Interaction 

1. Self-Trapping Transition 

Let us consider first the QMC results for the simplest case of the one-dimensional Holstein model. The variation 
of the ground state energy £^o(0) with the coupling constant A is shown in FIG. 01for a fixed phonon frequency of 
Oj — 1. The Eq{0) curve tends to the weak-coupling perturbation (WCP) result from below as A ^ 0, and to the 
strong-coupling perturbation (SCP) result as A ^ oo. Here we have also plotted the first and second terms of Eg. 1481 
which provide an indication of the potential energy (PE) and the kinetic energy (KE) of the system, respectively. 
Three separate regions can be clearly distinguished from FIG. 0J 

1. The large-polaron region at weak-coupling (A = to A w 1 in this case) is defined as the range of A for which 
Ef){0) is accurately described by WCP theory (section lV Ij() . At A = the energy is entirely kinetic and represents 
the bottom of the bare-electron band (with -Eo(O) = —zt). As A increases within this region, KE increases but 
remains large compared with PE (non- localized band- like states). 

2. The small-polaron region at strong-coupling (in this case A w 2.5 to A — > oo) has EqIO) accurately predicted by 
SCP theory (section IVAp . The PE is much greater than the KE (self-trapped states). 

3. The smooth transition region at intermediate coupling (A « 1 to A « 2.5 in this case) between the two regions 
above. As A increases in this region, a decrease in the PE and an increase in KE indicates the localization of 
the polaron. 
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FIG. 5: The variation of ground state energy £"0(0) (triangles), together with potential energy (diamonds) and kinetic energy 
(squares), with coupling A, for the one-dimensional Holstein model with a phonon frequency of lj — 0.5 (small symbols), 1.0 
(medium-size symbols) and 3.0 (large symbols). As u increases, the transition region shifts to higher A and becomes broader. 




FIG. 6: The number of phonons for the Holstein model at = 0.5, 1.0, 3.0 (small, medium-size, and large circles respectively). 
The dashed line shows the SCP result H57^ . The transition region boundaries are the same as those in FIG. 
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FIG. 7: The inverse effective mass mo/m* (where mo — /2ta^ is the bare-electron mass) for the Holstein model at lu — 
0.5,1.0,3.0 (small, medium-size, and large circles respectively). The dotted line shows the WCP result 18111 for the tight- 
binding Hamiltonian. 




FIG. 8: The logarithm of the effective mass for the Holstein model at w = 0.5, 1.0 and 3.0 (small, medium-size, and large circles 
respectively). The dashed line shows the SCP result ln(m*/mo) = 2X/ui. 
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FIG. 9: The isotope exponent for the Holstein model ai lo — 0.5, 1.0,3.0 (small, medium-size, and large circles respectively). 
The dotted line shows the SCP result. 

2. Variation with Phonon Frequency 

Now let us consider the way in which the properties of the Holstein polaron are affected by altering the value of 
the diniensionless phonon frequency a), as defined by Eq.|Sl This quantity (also known as the adiabatic ratio) is often 
used as a parameter in analytical approaches to the polaron problem: the adiabatic regime is defined as the case 
when a) < 1, and the anti-adiabatic regime as cZ; > 1. With this in mind, we present below the properties of the 
one-dimensional Holstein polaron for Cu — 0.5, uj — I, and iIi = 3. 

The ground state energies for all three values of w are presented together in FIG. against A. The PE tends to the 
same (w-independent) SCP result of Eo{0) = 2A as A ^ oo. Notice that the KE for intermediate and large values of 
A (mainly due to that of the lattice vibrations) decreases as u) increases. This affects the position of the start and end 
of the transition region. More precisely, as Co increases: 

1. The start of the smooth transition region moves to higher A. 

2. The transition region becomes broader (the start of the small-polaron region is also shifted to larger A). 

The values of A that mark the estimated start and end of the transition region are shown in TABLEIhII The definit ion 
of A, in Eq. is independent of to, and so naturally characterizes the three regions defined above. 

The variation of A^ph with A is presented in FIG. IHlfor each value of lu. Over all couplings, iVph decreases as the 
value of uj increases. Loosely speaking, it is simply "harder" to create phonons of higher frequency. The smooth 
transition from large to small polaron is again visible in the results for A^ph, with the edges of the transition region 
occurring at the same A as in the above results for Eo{0). 

The QMC results for the inverse effective mass mp / m* , and the isotope exponent on the effective mass , are 
presented against A in FIGS. [7H5| for each value of (D. We see from FIG. [7|that increasing co reduces the effective 
mass over the entire range of A. Both mo/m* and am* tend to the WCP solution as A becomes small, and to the 
SCP result as A becomes large. However, the results for ln(m*/mo) tend to the SCP solution at a slower rate than 
for the other observables. 

From all the results for the Holstein model above, we observe that as the value of ui increases, the curve for each 
observable (PE, A'ph, m*, and am») against A moves closer to the line representing the SCP result over the entire 
range of A. In other words, the SCP prediction becomes more applicable as cu increases. 
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A (end of large-polaron region) 


A (start of small-polaron region) 


0.5 


(1.1 ±0.1) 


(2.0 ±0.2) 


1 


(1.3 ±0.1) 


(2.7 ±0.3) 


3 


(2.1 ±0.1) 


(5.5 ±0.6) 



TABLE III: The boundaries of the transition region for the one-dimensional Holstein model at different values of the (dimen- 
sionless) phonon frequency lj. The estimates are based on the QMC results for the energies and A'^ph (observing the changing 
"trend" in plots of the deviation from the corresponding SCP and WCP result). 




FIG. 10: The ground state energy i5o(0) (triangles), potential energy (squares), and kinetic energy (diamonds) of the one- 
dimensional screened Frohlich model at = 1, versus A, for screening lengths R^c = 0, 1,3, oo (increasing size of symbols). 
The curves for £'o(0) (and PE) tends to the same strong coupling perturbation (SCP) result of Eo/t = — 2A (dashed line) as 
A — > 00. Note the crossing of the potential energy curves near A = 2. 

B. Screened Frohlich Interaction 

Here we consider the case of the one-dimensional screened Frohlich model, Eq. |H1 at a fixed dimensionless phonon 
frequency of u) = 1. We investigate the way the polaron properties depend on the range of the electron- phonon 
interaction by comparing the QMC results at four different values of the screening length (shown in FIG.|21): Rsc ~> 
(the short-range Holstein interaction), i?sc = !> ^sc = 3, and Rsc oo (the non-screened Frohlich interaction). 

The QMC results for the energy (i?o(0), PE, and KE) and the number of phonons in the polaron cloud N^i^ are 
presented against A in FIGS . [TUI and ITTI respectively. They are in excellent agreement with the WCP results at small 
A (not shown), and tend to the same (that is, 7- independent) SCP result as A ^ 00. One can see that, as the value 
of Rsc increases, the KE increases less rapidly with A, and so: 

1. The start of the transition region shifts to higher A. 

2. The transition region becomes broader (in A). 

This behavior is similar to that found for the Holstein model with increasing uj. In the present case, with the phonon 
frequency uj fixed, the above effect is due only to changing the shape of the electron-phonon interaction force /m(n). 
TABLE Hvl shows the values of A that mark the estimated start and end of the transition region for each value of i?sc- 
From the Afph results in FIG. IllL we notice that increasing the value of Rsc has the effect of increasing Afph at small 
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FIG. 11: The number of phonons in the polaron cloud A'^ph versus A for screening lengths Rsc ~ 0, 1, 3, co (increasing size of 
circles) at ij) = 1. The curves tend to the same SCP result of TVph — z\/oj (dashed line) as A ^ oo. 



Interaction Model 


A (start of transition) 


A(end of transition) 


Holstein (i?sc 0) 


(1.1 ±0.1) 


(2.0 ±0.2) 


Screened Frohlich (_Rsc ~ 1) 


(1.7 ±0.1) 


(3.8 ±0.4) 


Screened Frohlich (_Rsc ~ 3) 


(2.2 ±0.1) 


(4.7 ±0.4) 


Non-screened Frohlich {Rsc oo) 


(2.9 ±0.2) 


(6.3 ±0.5) 



TABLE IV: The boundaries of the transition region for the one- dimensional screened Frohlich interaction, at various screening 
lengths Rsc (measured in units of the lattice constant). The estimates are based on the QMC results for the energy and A'ph. 

A, and the opposite effect of decreasing Np^^ at large A. That is, the order that A^ph appears (with increasing i?sc) in 
the large polaron region is the reverse of the order in the small polaron region. 

The results for the effective mass are presented in FIGS. [T^ and [TSl and the isotope exponent results in FIG. 1141 
against A for the same four values of Rsc- We can see that for each value of Rsc the QMG results tend to the "model 
dependent" SGP result (of ln(m*/mo) = 2^X/Ld and am* — tA/cS, where the value of the mass enhancement factor 
7 is given in TABLE ^ for each i?sc)- This model-dependency is in contrast to the above results for Eo{0) and A^ph, 
which are 7-independent in the limit A — > 00. 

An important observation is evident from the plot of ln(m*/rno) against A shown in FIG. 1131 At intermediate and 
large couplings (that is, in the transition and small polaron regions), altering the value of -Rsc has a dramatic effect 
on the effective mass. For example, the non-screened Frohlich polaron is over 10"^ times "lighter" than the Holstein 
polaron at A = 4, and over 10^ times lighter at A = 5. 

The isotope exponent am- in FIG. 1141 shows a strong dependence on the range of electron-phonon interaction 
(as well as on uj and A) in the (physically most realistic) intermediate values of coupling. This is important, as 
experimental measurement (by Zhao et al.) of the exponent of the isotope exponent on the effective supercarrier 

mass along the Cu02 planes, jtt-q^, in the material La2-2;Sra;Cu04 shows a large value of a^^*^ = 1.9(2) in the deeply 

underdoped regime (x = 0.06), and a much smaller value of = 0.46(5) for optimal doping (x = 0.15)^"'^^. Both 
the magnitude and radius of the electron-phonon interaction should decrease with doping due to screening. These 
experimental results show that the former effect is more significant. 

Now let us consider the large-polaron region. The QMC results for to* and am* tend to the WCP results for all 
Rsc (not shown) as A ^ 0. As was the case for A^ph, we see from FIGS . [T^ and [T^ that the order that m* and am* 
appear in the large polaron region (with increasing Rgc) is the opposite to that in the small polaron region. As can 
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FIG. 12: The inverse effective mass mo/m* for screening lengtlis Rsc = 0, 1,3, oo (increasing size of circles) versus A at fixed 
(D = 1. For weak coupling (A < 1) the Holstein large-polaron has a slightly smaller m* than the long-range interactions. 




FIG. 13: The logarithm of the effective mass for screening lengths iisc = 0, 1, 3, oo (increasing size of circles) versus A at = 1. 
At intermediate and strong coupling, decreasing the value of Rsc dramatically increases the effective mass. The curves tend to 
the SCP result (dashed lines) at a slower rate than Eo{Q) and A/ph. 
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FIG. 14: The isotope exponent on the effective mass am* for screening lengths Rsc = 0, 1, 3, oo (increasing size of circles) versus 
A at a; = 1. 

be seen in FIG. IT^ the decrease of toq/to* with A is approximately hnear for the Holstein model, and approximately 
exponential for the other screening lengths. In fact, the effective mass of the Frohlich large-polaron is larger (up to 
approximately 10%) than the effective mass of the (short-range) Holstein large-polaron. 

It is apparent from the above results for the screened Frohlich model that as i?sc increases (from Holstein to 
Frohlich), the QMC results move, in general, closer to the SCP prediction over the entire range of A. That is, the 
SCP prediction becomes generally more applicable as the range of interaction increases (as well as with increasing 
uj). This fact allows us to define the intermediate region of the coupling strength and of the adiabatic ratio as a 
"small polaron" regime for any realistic-range e-ph interaction, that is the regime which is well described by the small 
polaron theory based on the Lang-Firsov transformed Hamiltonian averaged over phonons. 

C. Comparison with other approaches 

As a test of our method, we compare our results with those of other authors. Our main and original computations 
are for finite values of R; however, the majority of published work on lattice polarons relates to the Holstein {R = 0) 
interaction. FIG. ^1 compares our ground state energy for cu — 1.0, R — (compare FIG. 2)) with that obtained by 
other authors. To highlight the differences, the weak coupling approximation (|74|) is subtracted from the energy. The 
figure compares our results with the QMC data of Hohenadler et al.'^^ (based on a bosonic path integral, evaluated 
at inverse temperature /3 = 10) and exact diagonalization results for small clusters££. Different computations yield 
similar values of the ground state energy, but our QMC energies are closer to the the exact diagonalization results; 
our calculations have been carried out at a lower temperature, /3 — 25. FIG. \W\ shows good agreement between our 
results for the effective mass (compare FIG.|Hll variational calculations of Bonca et alM- (Note that the A in that 
work corresponds to \/2ujt\ in our notation.) We also reinforce earlier conclusions on the dependence of effective 
mass with interaction range^Si^ by interpolation between the Holstein and Frohlich limits. 
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FIG. 16: ID Holstein polaron effective mass for ui = 1.0, R = (crosses) compared with variational results (solid line, Ref.^^)- 
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VII. CONCLUSIONS 

The general aim of this work was to investigate the way in which the range of electron-phonon interaction governs 
the physical properties of the (single) lattice polaron. The understanding of this is of considerable current interest 
because of the increasing amount of experimental evidence suggesting that polarons are present in high-temperature 
superconducting and colossal magnetoresistance materials. 

Perturbation approaches are used in the limits of strong and weak electron-phonon coupling strength A. However, 
in general these do not provide an acceptable description in the (physically most realistic) intermediate coupling range 
A « 1 . We have performed an extensive Monte Carlo study of the ground state properties for the screened Frohlich 
polaron, equation (jS]), in one dimension, over a wide range of coupling. 

We have used path-integral quantum Monte Carlo, in which the phonon degrees of freedom are analytically inte- 
grated out, leaving only the electron coordinates to be simulated. The use of a path integral with twisted boundary 
conditions allowed us to extract dynamic properties directly from the simulations. There were no (systematic) errors 
due to finite size or finite time-step. 

The properties measured were the ground state energy Eq{0), the number of phonons in the polaron cloud iVph, 
the effective mass m*, and the isotope exponent on the effective mass am*- The QMC results were always found to 
tend to the weak-couphng perturbation (WCP) predictions for A — > 0, and to the strong-couphng perturbation (SCP) 
predictions for A — *■ oo. 

The screened Frohlich polaron was studied for various values of the screening length Rac (which essentially controls 
the range of the electron-phonon interaction): R^c ^ (on-site, Holstein interaction), i?sc = 1, -Rsc = 3, and 
Rsc oo (long-range, non-screened Frohlich interaction). For each value of Rsc, we determined the variation of the 
above observables with A, at a fixed phonon frequency of w = 1. The main findings are summarized below: 

1. We observe the presence of a self-trapping transition for all values of Rsc- In each case, the following three 
regions are identified: 

(a) The large-polaron region at weak coupling, in which the behavior of the system is accurately described by 
WCP theory. This region is characterized by delocalized, band-electron-like states. 

(b) The small-polaron region at strong coupling, in which the behavior is accurately described by SCP theory. 
This region is characterized by localized ("self-trapped") polaronic states. 

(c) The transition region between the two, at intermediate coupling. We observe a smooth crossover from large 
to small polaron in all the observables measured. 

2. The transition region boundaries depend on the range of interaction. As the value of i?sc increases 

(a) The start of the transition region (the point at which it becomes energetically favorable for localized states 
to exist) shifts to higher A. 

(b) The transition region becomes broader (the start of the small polaron region also shifts to higher A). The 
small-polaron region starts when the kinetic energy is much smaller than the potential energy. 

(c) The values of the observables (PE, A^ph, rn*, and am') generally move closer to the corresponding SCP 
result over the entire range of A. 

3. In the large polaron region, the effective mass for long-range electron-phonon interaction (i?sc > 1) is found to 
be up to approximately 10% larger than that for the Holstein interaction {Rsc — *■ 0). 

4. We observe large variations in the isotope exponent on the effective mass am* in the (physically most realistic) 
intermediate coupling regime (with changing i?sc and A, as well as lu). This is encouraging, as experimental 
observation shows large variations in the isotope exponent with the level of doping in high-Tc materials. 

5. Reducing the range of electron-phonon interaction dramatically increases the effective mass for intermediate 
and large values of coupling (that is, in the transition and small-polaron regions). In comparison, Eq{0) and 
TVph are only slightly affected by altering Rgc- 

6. We also study the dependence on phonon frequency oj of the Holstein polaron {Rsc ^0). As w increases 

(a) The transition region shifts to higher A. 

(b) The transition region becomes broader. 

(c) The observables (PE, Aph, m* , and am*) move towards the corresponding SCP result over the entire range 
of A. 
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7. Increasing interaction range has a qualitatively similar effect to increasing phonon frequency; compare for 
example FIGS. El and CI or FIGS. 1^ and [Tl 

One can also calculate the isotope effect on the whole polaron band dispersion applying the continuous-time quantum 
Monte-Carlo algorithm^^. To deal with the electron spectral function and high-energy excitations, involving phonon 
shake-off, measured in ARPES^SiiMi, the QMC algorithm has to include the off-diagonal paths, which remains a 
challenging but solvable problem of our QMC simulations in the site representation. Also other methods as the 
momentum based QMGiS, the numerical diagonalization of vibrating clusters or the renormalisation group are able 
to calculate the spectral function. The method used here relies on a phonon gap, the errors being exponentially small 
in flhuj, and is therefore more accurate for intermediate and large values of to. While lower lu can be also simulated 
without major difficulty but with increased inverse temperature and therefore increased CPU time, the parameters 
used cover the most physically relevant rangA 
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